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ABSTRACT 


The  purpose  of  the  present  investigation  is  the  deter¬ 
mination  of  the  properties  of  an  intense  10.6-micron  laser 
beam  pi  opagating  through  the  open  atmosphere  in  the  pres¬ 
ence  of  wind  or  slewing,  or  both.  It  is  shown  that  the  Max¬ 
well  equations  can,  in  this  problem,  be  reduced  to  the  study 
of  the  scalar  wave  equation  with  a  varying  index  of  refrac¬ 
tion.  The  index  of  refraction  is  related  to  the  atmospheric 
density;  therefore,  the  density  changes  in  the  air  due  to 
beam  absorption  are  related  to  the  absorption  coefficient  of 
the  air  and  to  the  intensity  of  the  beam,  using  the  linearized 
hydrodynamic  equations.  A  detailed  discussion  of  the  mech¬ 
anisms  of  photon  absorption  by  the  constituents  of  the  air  is 
presented.  The  resultant  equation  for  the  scalar  wave  is  a 
nonlinear  partial  differential  integral  equation  which  is 
solved  numerically.  The  algorithm  used  for  the  computer 
code  is  discussed,  together  with  criteria  that  have  been  de¬ 
termined  to  be  useful  in  assessing  the  accuracy  and  relia¬ 
bility  of  the  numerical  results.  The  solutions  of  several 
different  problems  are  presented  and  discussed.  In  partic¬ 
ular,  it  is  found  that  (a)  beam  quality  is  degraded  for  water 
vapor  pressures  at  or  near  sea  level,  and  (b)  beam  slewing 
reduces  the  detrimental  effect  of  water  vapor. 
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PROPAGATION  OF  HIGH-ENERGY  10.6-MICRON  LASER  BEAMS 
THROUGH  THE  ATMOSPHERE 


I.  INTRODUCTION 

The  development  of  the  gas  dynamic  laser  has  turned  the  prospect  of  a  cw  high- 
power  laser  into  a  reality.  Many  applications  of  these  devices  involve  propagation  of  the 
beam  through  the  atmosphere,  and  thus  an  understanding  of  the  various  aspects  of  beam- 
atmosphere  interactions  is  required.  This  report  is  devoted  to  a  study  of  one  of  the  non¬ 
linear  aspects  of  the  propagation  of  laser  radiation  through  air,  namely,  the  change  due 
to  absorption  of  beam  energy  in  the  index  of  refraction  of  the  air,  which  in  turn  modifies 
the  propagation  of  the  beam.  Beam  energy  absorption  by  an  otherwise  motionless  atmos¬ 
phere  induces  temperature  changes  in  the  medium,  which  in  turn  cause  a  defocusing  (or 
focusing)  of  the  beam  that  changes  with  time.  A  number  of  studies  have  been  made  of 
this  phenomenon,  which  is  often  referred  to  as  thermal  blooming  (1-7). 

If  a  steady  wind  is  present  with  a  nonzero  component  transverse  to  the  beam  axis, 
the  heated  air  is  swept  out  of  the  beam  and  a  steady-state  ultimately  evolves  (8).  Beam 
slewing  under  suitable  conditions  likewise  gives  rise  to  a  steady  state.  The  formulation 
and  description  of  these  steady  states  is  the  problem  this  report  addresses.  The  exist¬ 
ence  and  stability  of  such  states  is  assumed.  Theoretical  and  experimental  studies  of 
cases  such  as  these  have  been  made  (9). 

Section  II  of  this  report  presents  the  formulation  of  the  physical  assumptions  of  the 
problem.  Section  HI  is  a  discussion  of  the  techniques  used  in  a  numerical  solution  of  the 
equations  derived  in  Section  n,  while  Section  IV  is  a  discussion  of  the  results  of  several 
numerical  solutions,  and  a  presentation  of  the  basic  limitations  of  the  algorithm  used  to 
obtain  the  numerical  solutions. 


II.  FORMULATION  OF  THE  PROBLEM 
Outline 

An  equation  which  contains  the  index  of  refraction  of  the  air  is  developed  to  describe 
the  propagation  of  the  laser  beam.  The  index  and  the  density  of  the  air  are  related  by 
the  Lorentz-Lorenz  law.  The  air  density  is  in  turn  related  to  the  heat  sources  in  the  air 
that  arise  from  beam  absorption  through  the  linearized  hydrodynamic  equations.  Finally, 
a  detailed  discussion  of  the  mechanism  of  absorption  of  10.6-micron  radiation  along  the 
lines  developed  by  Wood,  Camac,  and  Gerry  (10)  is  included. 


Propagation  Equation 

The  wave  equation  for  a  harmonically  varying  electric  field  in  a  medium  of  varying 
dielectric  coefficient  is 


V2e  +  k  2«E  =  W  •  E 


1 
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with  k  =  u  'c ,  where  u  is  the  angular  frequency,  and  e  is  the  dielectric  coefficient.  Us¬ 
ing  the  equation  of  charge  conservation,  the  divergence  term  on  the  right-hand  side  is 
given  by 


V  •  E  =  -  <T 1  E  •  Ve  . 

If  the  properties  of  the  medium  are  such  that  variations  of  the  dielectric  coefficient  are 
small  in  a  wavelength,  the  right-hand  term  may  be  neglected  in  comparison  with  the  k2tE 
term  in  the  first  equation,  yielding  simply 


VJE  +  k  2«E  =  0  . 


This  approximation  also  implies  that  the  polarization  of  the  field  is  not  changed  as  it 
propagates  so  that  only  linearly  polarized  amplitudes  need  be  considered.  Further,  since 
the  principal  concern  is  with  a  beam  which  in  the  first  approximation  is  simply  propa¬ 
gating  along  the  z  axis,  a  solution  of  the  form  E  =  e<p(r)  exp(ik  sJT^  z)  is  anticipated, 
where  <p  is  assumed  to  be  a  slowly  varying  function  of  z ,  c0  is  the  dielectric  constant  of 
the  medium  prior  to  its  interaction  with  the  beam,  and  e  is  a  unit  vector  transverse  to 
the  beam  axis.  £  and  e0  each  have  an  imaginary  part  that  describes  beam  absorption; 
because  density  changes  will  be  small,  =  c0(i)  to  an  excellent  approximation.  There¬ 
fore,  the  equation  for  v  is 

—  +  2ik  V  2  (p  +  k2(£(r>  -  e0(r))<P  =  0 
3z2  <)z  1 

where  32/3x2  t  3  V  3y  2 .  in  accord  with  the  expectation  that  <j>  be  a  slowly  varying 
function  of  z,  the  second-derivative  term  is  dropped  as  being  small  compared  with  the 
first-derivative  term.  The  equation  that  describes  the  laser  beam  then  is 

2ik  — +  t  k2(£<r>-€0<r>)v  =  0  .  (1) 


Equation  (1),  being  essentially  a  two-dimensional  Schr Winger  equation  with  z  playing  the 
role  of  time,  has  the  quantity 


i= 


If 


<P*(  x,  y,  z)  <p(x,y,z)  (lx  dy 

ant 


(2) 


as  a  constant  of  the  motion. 

The  energy  flux  is  described  by  the  time-averaged  Poynting  vector.  With  the  ap¬ 
proximations  described  above  included  in  the  standard  expression  for  the  Poynting  vector, 
the  energy  flux  becomes 


P  =  £  e—  ¥**  ?z 


(3) 


where  a  =  2€0(i)  .  The  total  power  P  passing  through  a  plane  z  =  constant  is  given  by 


P(z)  =  e_<XZ  JJ  dx  dy  <p*?> 


t»  const  ant 


(4) 
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e z  is  unit  vector  along  z  axis.  Equation  (1),  supplemented  with  a  relationship  between  the 
dielectric  constant  and  the  function  <p,  is  the  equation  to  be  solved. 


Hydrodynamics 

The  atmosphere  will  be  treated  as  a  perfect  fluid;  viscous  and  thermal  conduction 
effects  will  be  ignored  since  their  consequences  can  be  shown  to  be  negligible  in  the 
cases  that  will  be  considered  in  this  paper.  The  beam  will  be  regarded  as  an  external 
heat  source  described  by  an  energy  deposition  function  Q ,  representing  the  energy  de¬ 
posited  per  gram  per  second  at  each  point  in  the  fluid.  A  complete  description  of  the 
hydrodynamic  system  is  thus  provided  by  the  boundary  conditions,  the  equation  of  state, 
an  expression  for  the  internal  energy  of  the  fluid,  and  the  hydrodynamic  equations  of 
motion: 


0  . 


0  . 


pQ 

Cyf  , 


(cp-cv)T  . 

The  quantities  p,  y,  p,  g,  T,  cv,  and  cp  represent  the  density,  fluid  velocity,  pressure, 
internal  energy  per  unit  mass,  temperature,  specific  heat  at  constant  volume,  and  spe¬ 
cific  heat  at  constant  pressure,  respectively.  The  operator  d/dt  stands  for  the  sub¬ 
stantial  derivative. 

During  the  course  of  fluid  motion  through  the  laser  beam,  the  temperature  ch£U  ge 
of  a  fluid  element  will,  in  most  circumstances,  be  small  so  that  the  specific  heats  may 
be  regarded  as  constant  throughout  the  flow.  Under  these  conditions,  the  first,  third, 
fourth,  and  fifth  of  these  equations  may  be  combined  to  yield  an  expression  involving 
pressure  and  density  alone: 


af  +  ^‘v  = 

dg  n 

p  __  +  pV  •  y  = 

g  = 

P  _ 
P  ' 


dp  p  do 


where  y  -  cp/cv.  Further,  because  the  changes  in  the  pressure,  density,  and  velocity 
are  expected  to  be  small,  the  linearized  hydrodynamic  equations  can  be  regarded  as  a 
good  description  of  the  physical  system.  Let  the  pressure  p0 ,  density  p0 ,  and  velocity 
v0  of  the  fluid  satisfy  the  unperturbed  exact  hydrodynamic  equations 


d£0 

dt  + 


=  0 


dv0  - 

Po  "dF  +  =  0  • 

dP0  Po  dPo  Q 
"dt  "  y  ~pl  "dt" 
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The  total  pressure,  density,  and  velocity  are  taken  to  be 

P  =  P0  +  Pi  • 

P  -  Po  f  Pi  • 


The  subscripted  quantities  p1(  Pi,  and  v,  are  taken  to  be  of  the  first  order,  and  higher 
powers  of  these  terms  are  deemed  negligible.  Particularizing  to  the  case  of  a  beam  in 
the  presence  of  a  steady  uniform  wind  ( v0  =.-  constant),  the  first-order  quantities  are 
taken  to  satisfy  the  equations 


dp.  _ 

-fr*  Po?-v,  =  0 

(5) 

dv,  _ 

Po  TiF  +  ^  i  =  0 

(6) 

dp,  Po  dPl 

dt  ‘  yJ0  -dT-  poQ 

(7) 

where,  now,  the  operator  d/dt  is  taken  to  be 

d  3  3 

dt  =  It  +  v«  '  V  ' 

(8) 

The  term  yp0/p0  is  cs2,  where  cs  is  the  speed  of  sound  in  the 
tions  caused  by  the  beam. 

fluid  prior  to  the  perturba- 

Since  the  steady-state  case  is  the  solution  of  interest,  the  hydrodynamic  quantities 
at  a  fixed  point  in  space  are  taken  to  be  independent  of  time.  Therefore  all  partial  time 
derivatives  of  these  quantities  vanish,  and  the  equations  of  motion  reduce  to 

-*  -♦ 

V0  •  VP 1  +  Pov  •  V,  =  o  , 

(9) 

PovO  •  Vv,  +  Vp,  r  0  , 

(10) 

V0  'V  (Pi-Cs^i)  =  (7-  !)  PoQ  • 

(11) 

The  last  equation  may  be  integrated  to  give 

(r-i)  f° 

v„c  2  J 

V0'-J  -® 


Pi  = 


ds  p0Q(  r  +  sv0)  + 


Pi 


(12) 


where  v0  is  a  unit  vector  along  the  direction  of  the  wind  velocity  and  the  integral  is 
along  the  path  of  motion  of  a  fluid  element.  The  path  is  parameterized  by  the  path  length 
s  which  is  taken  to  be  zero  when  the  fluid  element  is  at  the  point  r .  By  elimination  of 
the  density  and  velocity  in  favor  of  the  pressure,  a  differential  equation  for  the  pressure 
change  alone  is  easily  obtained: 
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Pi  = - —  vo  •  J  f  o'? 


(13) 


An  estimate  of  the  contribution  of  the  pressure  term  to  the  density  can  be  obtained  by 
ignoring  the  pressure  gradient  along  the  beam  axis  and  the  gradient  perpendicular  to  the 
wind  velocity;  then 


p  j(  r)  =  - 


(y-  i)  v0 


4 


J  ds  o0  Q 


9  (r>svc)  . 


(14) 


Restricting  wind  speeds  to  values  much  smaller  than  the  speed  of  sound,  the  pressure 
term  may  clearly  be  ignored,  and  the  density  change  of  air  at  any  point  r  is  given  by 


(y-  1)  f 

Pi(  r)  - - —  I  <ls  p0Q(r  +  SV0)  . 

v«c. 


(15) 


Equation  (15)  shows  clearly  that  a  zero  wind  speed  is  not  tenable  in  this  treatment;  a 
steady  state  cannot  be  achieved  without  the  wind  as  a  mechanism  for  removing  the  heated 
air  from  the  beam.  Furthermore,  Eq.  (15)  can  be  used  as  a  means  for  estimating  the 
limits  that  must  be  placed  upon  the  beam  intensity,  wind  velocity,  and  absorption  coef¬ 
ficient  for  the  linearized  treatment  used  here  to  retain  its  validity.  If  beam  absorption 
leads  instantaneously  to  heating,  then  the  quantity  p0Q  becomes  ai ,  where  I  is  the  beam 
intensity.  The  integral  may  be  approximated  by  2oIa  where  I  is  an  average  power  dis¬ 
tribution  along  a  line  across  the  laser  face  through  the  center,  and  a  is  the  beam  radius. 
Then,  if  the  density  change  must  be  less  than  p ,  to  preserve  the  validity  of  the  lineariza¬ 
tion  of  the  hydrodynamic  equations,  1,  «,  and  v0  must  together  satisfy  the  inequality 


a  In 


filcs 

2(y-l) 


(16) 


In  this  formulation  of  the  propagation  of  a  laser  beam  in  the  presence  of  a  steady 
wind,  there  exists  a  reflection  symmetry  in  the  fluid  variables  and  the  beam  intensity 
through  a  plane  encompassing  the  beam  axis  and  a  line  parallel  to  the  wind  velocity. 

This  symmetry  has  considerable  utility  in  the  numerical  solution  of  the  problem. 

Important  applications  may  be  envisaged  wherein  the  laser  beam  is  required  to 
rotate  about  an  axis  perpendicular  to  the  beam  axis;  such  motion,  referred  to  henceforth 
as  beam  slewing,  can  be  brought  within  the  framework  of  a  steady-state  calculation  sim¬ 
ilar  to  the  discussion  above.  Beam  slewing  in  an  otherwise  stationary  atmosphere  will 
not  lead  to  a  steady  state,  in  part  because  there  is  no  mechanism  for  removing  the  heated 
air  in  the  vicinity  of  the  laser  face,  even  though  the  slewing  does  so  further  downrange. 
Therefore,  the  initial  conditions  will  be  functions  of  time.  A  steady  state  can  be  achieved 
by  the  introduction  of  a  mechanism  for  removing  the  heated  air  from  the  vicinity  of  the 
face  of  the  laser;  such  a  mechanism  may  be  modeled  in  a  wide  variety  of  ways,  three  of 
which  will  be  briefly  discussed  here. 

In  a  coordinate  system  rotating  with  the  laser,  the  hydrodynamic  equations  assume 
the  form 


(&•  v4 


-4  ■ 


(17) 
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<  v  '  vjv  =  -  i\7p  -  2J2  x  v  -  fix  (fix  r)  (18) 

(s,v'7):7(rv'v}t^1^  d9) 

where  v  is  the  fluid  velocity  measured  relative  to  the  rotating  reference  frame,  p  and  H 
represent  the  pressure  and  density,  respectively,  while  n  is  the  angular  slewing  rate, 
which  is  orthogonal  to  the  beam  direction.  These  equations,  when  linearized,  assume  the 
forms,  for  the  steady  state, 


(v0  •  7 )p,  =  -  p0V  ■  V, 

(20) 

(v  ■  V)v,  +  0  x  V,  =  -  ~  3p, 

P  0 

(21) 

<vo  •?)  <Pi-csVi)  -  (y-t/PoQ 

(22) 

while  the  zeroth-order  quantities  satisfy  equations  of  the  form  of  Eqs.  (17)- (19)  with  the 
heating  term  dropped.  An  exact  solution  of  the  (nonlinear)  zeroth-order  equations  is 
provided  by 

V0  =  Vf0  -  0  X  r  (23) 

where  *0  is  a  velocity  vector  parallel  to  Q.  For  nonvanishing  w0,  this  solution  corre¬ 
sponds  to  a  wind  along  the  slewing  axis,  and  this  wind  serves  as  a  mechanism  for  remov¬ 
ing  the  heated  air  from  the  laser  face.  In  the  rotating  coordinate  system,  a  fluid  particle 
is  moving  in  a  spiral  path  about  the  slewing  axis,  and  its  trajectory  is  given  by 


r(  s) 


Os  .  Qs\. 

yo  cosT7+  z«  s,nv7jJ 


I  ■  Os  Os\,  .... 

\  y°  sin  ^7  *  z°  cos  ^jk  (24) 

while  the  magnitude  of  the  fluid  velocity  is  (w02  ^  02  z2)17  2  and  is  a  constant  for  a  given 
fluid  element.  Under  these  conditions,  Eq.  (22)  may  je  integrated  to  give 


t> i<  K  s)) 


(y-i) 


p0Q(r(s‘)) 


P  i(  K  s>) 


(25) 


The  point  r0  may  be  chosen  to  represent  a  point  in  space  where  the  beam  has  not 
arrived,  so  that  the  density  and  pressure  changes  there  vanish;  hence  the  factor  in 
brackets  in  Eq.  (25)  may  be  set  equal  to  zero.  An  expression  for  the  pressure  change 
alone  may  be  derived  in  much  the  same  manner  as  before.  Again,  the  contributions  of 
the  pressure  term  can  be  shown  to  be  negligible  provided  that  v0  •  <  cs;  however,  with 
slewing  present,  v0  can  become  quite  large  for  long  distances  down  the  beam.  There¬ 
fore,  the  linearization  of  the  hydrodynamic  equations  will  be  valid  only  for  restricted 
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distances  zM  downbeam  for  which  zM  «  cs/n.  Within  these  limitations,  the  solution  for 
the  density  change  is  given  by 


s 

Pi(  r)  =  -  — — —  f  ds'  p00(r(s'))  , 
_  0  J{\ 


v0  =  (w02  f  n2z2)  .  (27) 

On  the  basis  of  the  above  solution,  three  models  can  be  described. 

1.  The  velocity  component  w0  is  nonvanishing.  As  previously  mentioned,  a  wind  of 
this  type  serves  as  a  device  for  removing  the  heated  air  from  the  laser  face,  allowing  a 
steady  state  to  be  achieved,  and  permits  an  arbitrary  number  of  rotations  of  the  laser 
about  its  slewing  axis.  This  solution  lias  the  disadvantage  that  the  symmetry  described 
at  the  end  of  the  formulation  of  the  steady-wind  solution  is  lost,  and  therefore  doubles 
the  computational  effort  in  a  numerical  solution  of  the  problem. 

2.  Set  the  velocity  component  w0  =  0,  and  do  not  allow  the  axis  of  revolution  to  go 
through  the  laser  face,  but  instead  let  it  be  at  a  distance  d  behind  it.  Then  the  rotation 
itself  removes  the  heated  air  from  the  laser  face;  all  that  is  required  is  that  the  steady 
state  set  in  prior  to  one  complete  rotation  of  the  laser  about  the  slewing  axis.  If  a  dis¬ 
tance  R  down  the  beam  axis  is  measured  in  terms  of  z,  the  distance  from  the  laser  face, 
then 

V0  =  fi(D+  z) 

=  v;  t  nz  .  (28) 

This  method  has  the  advantage,  from  the  standpoint  of  numerical  computation,  of  pre¬ 
serving  the  symmetry  that  reduces  computation  times. 

3.  Again,  set  the  velocity  component  w0  =  0;  in  the  rotating  reference  frame, 
rigidly  attached  to  the  laser,  set  up  a  fan  that  blows  air  at  a  velocity  v'  at  right  angles 
to  the  beam  axis  and  the  slewing  axis,  in  addition  to  the  slewing  velocity,  for  a  fixed  dis¬ 
tance  D  down  the  beam  axis.  Then  v0  in  Eq.  (27)  is  replaced  by 


Vg  +  Qz  ,  0  <  z  <  D 


which,  in  for  a,  is  similar  to  the  second  model  for  slewing.  This  model  also  preserves 
the  desired  symmetry. 

The  model  used  in  calculations  here  will  be  the  second  model.  All  three  models 
admit  a  new  parameter,  in  addition  to  the  slewing  rate  ft;  solutions  will  therefore  be 
dependent  upon  the  value  of  the  parameter.  This  parameter  dependence  is  expected  to 
be  discernible,  but  not  of  vital  importance. 
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Absorption  of  10.6-Micron  Radiation  By  Air;  Kinetic  Cooling 

The  constituents  of  air  that  are  the  principal  absorbers  of  10.6-micron  radiation  are 
water  vapor  and  carbon  dioxide.  Pressure  broadened  levels  of  water  vapor  make  it  an 
effective  absorber  of  infrared  radiation  throughout  a  broad  band  that  includes  10.6 
microns,  while  carbon  dioxide  absorbs  principally  through  the  inverse  laser  transition. 
These  processes  have  been  briefly  discussed  by  Wood,  Camac,  and  Gerry  (10)  and  have 
been  applied  in  a  study  of  C02  laser  propagation  in  geometrical  optics  by  Wallace  and 
Camac  (11).  Because  of  its  importance  to  the  character  of  the  results,  a  detailed  dis¬ 
cussion  of  the  absorption  mechanism  is  presented  here. 

Water  molecules,  upon  absorption  of  one  photon  of  10.6-micron  radiation,  are  ex¬ 
cited  to  the  (010)  vibrational  state.  The  stored  energy  is  released  to  the  translational 
modes  by  collisional  deexcitation  of  this  state;  the  principal  molecules  in  the  deactiva¬ 
tion  collisions  are  oxygen  and  water  vapor.  The  deactivation  time  is  so  short  compared 
to  all  hydrodynamic  processes  in  the  propagation  problem  that  it  may  be  regarded  as 
zero.  Hence,  the  energy  deposition  rate  may  be  written  as 


/>0Q(r,t)  =  aH2oI<r>t)  4  (30) 

where  p0Q'  is  the  energy  deposition  rate  due  to  all  other  constituents  of  air. 

Figure  1  is  an  energy  level  diagram  showing  the  principal  lower  vibrational  levels  of 
C02,  N2,  and  H20.  The  excited  levels  of  the  various  molecules  are  populated  according 
to  Boltzmann  statistics;  therefore  C02  molecules  in  the  (100)  state  are  always  present 
and  bean  photons  may  be  absorbed,  leaving  the  C02  molecules  in  the  (001)  state. 

The  reaction 


CO2(001)  +  N2  ;=±  C02(000)  +  N*  +  18  cm'1  (31) 

is  very  rapid  due  to  resonance  between  the  levels  of  C02  and  N*,  the  characteristic  time 
being  of  the  order  of  10~7  sec.  The  system  of  carbon  dioxide  molecules  is  thus  dis¬ 
placed  from  a  thermal  equilibrium  distribution  with  an  excess  of  molecules  in  the  ground 
state  and  the  (001)  state,  and  a  deficiency  in  the  (100)  level.  The  latter  state  will  be  re¬ 
populated  principally  by  collisions  of  the  type 

C02(010)  +  CO2(010)  ;=±  CO2(100)  +  C02(000)  -  54  cm'1  .  (32) 

This  reaction  proceeds  faster  than  all  others  that  tend  toward  the  same  final  state  be¬ 
cause  of  resonance.  The  energy  removed  in  this  reaction  comes  from  the  translational 
modes.  Reaction  (32)  tends  to  deplete  the  supply  of  (010)  levels,  but  these  are  replen¬ 
ished  through  the  process 

C02(000)  +  C02  (000)  ;=±  C02  (000)  +  CO2(010)  -  667  cm'1  .  (33) 

The  absorbed  energy  in  Reaction  (33)  once  again  comes  from  translation.  Two  reactions 
of  type  (33)  must  occur  for  every  one  of  the  type  indicated  by  Reaction  (32)  to  maintain 
the  C02  in  thermal  equilibrium.  The  removal  of  energy  from  the  translational  modes  by 
Reactions  (32)  and  (33)  cools  the  C02  molecular  system,  and,  concomitantly,  the  air. 
Because  the  reactions  subsequent  to  the  photon  absorption  that  are  described  above  occur 
so  rapidly  as  to  be  almost  instantaneous,  the  rate  of  energy  removal  from  the  transla¬ 
tional  modes  is  governed  by  the  rate  at  which  photon  absorption  occurs.  Therefore  the 
contribution  of  the  C02  molecules  to  p0Q'  is  given  by 
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Fig.  1  -  The  principal  lower  vibrational  energy  levels 
of  the  major  constituents  of  air— CO 2,  N2,  02,  and 
H20— are  shown,  with  the  10.6-micron  laser  transition 
explicitly  marked.  The  10.0-micron  photon  energy  is 
designated  as  . 


1 00 


=  “C021(r't) 


(34) 


Reaction  (31)  is  populating  the  N2  above  the  equilibrium  level;  simultaneous  with  the 
CO 2  collisions,  N*  is  undergoing  collisions  of  the  form 


+  02 


N,  +  Oj  +  £00,  -  18  cm'1 


and 


(35) 


001 


n£  +  h2o 


N2  +  H20  +  e 


-  18  cm' 1 
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which  tend  to  depopulate  the  excited  level  N*.  In  addition  to  Reactions  (35),  depopulation 
may  go  via  excited  states  of  O2  and  H2O,  and  Reactions  (35)  will  be  understood  to  sym¬ 
bolize  both.  When  water  vapor  is  present  in  any  significant  amount,  deactivation  of  nitro¬ 
gen  proceeds  primarily  along  the  second  route;  since  the  water  vapor  content  of  air  is  a 
variable  quantity,  both  reactions  must  be  included. 

Let  SN*  represent  the  excess  number  of  Nf  molecules  at  time  t  over  the  equilibrium 
number  at  some  point  r  in  the  beam.  Its  time  rate  of  change  is  given  by 


_d 

dt 


SN*  = 


aco2 

KUy 


I(r.t) 


I^NjOj-NjOj^Oj  +  ^N*Hj0-.N2H20NH20  j  (36) 

where  the  quantity  p  is  the  reaction  constant  for  the  reactions  indicated  by  the  subscripts 
(12),  and  No2  is  the  number  density  of  02  molecules,  and  NH20  the  number  density  of 
water  molecules.  In  terms  of  a  reaction  time  r,  the  above  equation  may  be  rewritten  as 


I  ( r, t ) 


T 


(37) 


where 


l  _  i  P(Oj)  j  p(H20) 

r  tn«o2  p  tn2h2o  P  ’  (38) 

p(02)  and  p(H20)  are  the  partial  pressures  of  O2  and  H20,  respectively,  P  is  the  total 
pressure,  and  tn$o2  and  rN.Hj0  are  the  reaction  times  for  Reactions  (35).  The  deriva¬ 
tive  is  the  substantial  time  derivative  and,  for  the  steady-state  case,  the  solution  to  the 
above  equation  becomes 


SN*(r) 


“CO, 


S 


v0fiw  -co 


J  ds  I(rt  svQ  )  eV° 


(39) 


The  number  of  energy-releasing  transitions  per  second  is  SN */r  and  the  rate  of 
energy  release  is 


001  C°2  f ,  w  *  ^  ^ 

_  __  J  Js  SV0) 


■ftu  v0  - 

y  “03 


(40) 


The  total  energy  deposition  rate  is  obtained  by  combining  Eqs.  (30),  (34),  and  (40) 

„  ,  €ioo  \ 

P°Q  =  ( ““2° '  ^7  ^c°2  J I(r) 


€  ®co  r° 
*001  I 

Ifoy  V0T  L 


ds  I(r+svQ)  e 


(41) 
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This  last  result,  taken  together  with  the  steady-state  solutions  of  the  linearized  hydro 
dynamic  equations,  yields,  after  evaluation  of  some  integrals, 


p,(r)  = 


Pifr) 


I 

1  where 


<7-  1) 


2 

s 


a 


I  (  r  +  s  v0  ) 


a 


co. 


+  a 


HjO 


and 


S  -  £°01  gC°2 

•h&iy  a 

In  addition,  the  fact  that 

€00i  -  €too  1 

has  been  incorporated  into  Eq.  (42). 


(42) 


(43) 


(44) 


(45) 


The  quantities  aCOj  and  a„2o  for  air  are  related  (13,  14)  to  the  partial  pressures  of 
C02  and  water  vapor,  temperature,  and  total  pressure  by 


a 


co. 


1.44  XlO'3 


970 

10  T  cm"1 


(46) 


and 


®H20  =  4.32  Xl0'up(H2O)!P+  193p(i!20>)  cm"1  .  (47) 

The  absorption  coefficient  of  water  vapor  in  air  used  here  has  been  determined  empiri¬ 
cally  only  in  the  temperature  range  of  23  °C  to  26 °C;  therefore,  Eq.  (47)  has  a  limited 
temperature  range  of  validity  that  must  be  borne  in  mind. 

If  the  deactivation  time  r  of  nitrogen  were  zero,  then  there  would  be  no  delay  be¬ 
tween  the  time  of  photon  absorption  and  the  appearance  of  the  photon  energy  as  heat. 

This  is  manifest  in  Eq.  (42)  where  the  term  in  the  integrand  that  involves  r  would  not  be 
present,  and  the  density  change  px  would  always  be  negative,  i.e.,  the  air  always  expands 
in  the  presence  of  the  beam.  Since  r  in  general  does  not  vanish,  varying  amounts  of  heat¬ 
ing  and/or  cooling  may  take  place.  For  the  case  of  no  water  vapor  in  the  air  at  all, 
o(h2o  =  0  and  $  =  2.441,  while  7  assumes  its  maximum  value.  Thus,  for  large  velocities, 
i.e.,  for  v0  such  that  »0«a/r  throughout  the  beam,  the  factor  (l  -  >>)  in  the  integrand  multi¬ 
plying  the  intensity  becomes  -1.44,  so  that  px  >  0  and  the  beam  absorption  causes  cool¬ 
ing.  For  lesser  wind  speeds  and  increasing  water  vapor  content,  the  physical  situation 
will  be  intermediate  to  these  two  extremes,  but  it  should  be  noted  that  there  can  be  no 
net  cooling  when  aH  0  >  1.44  aCOj. 
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The  density  change,  as  expressed  in  Eq.  (47)  (with  the  pressure  term  neglected),  is 
related  to  the  wave  function  <p  by  Eq.  (4): 


P,(r)  =  -  o-aj 

voc. 


J  ds  (l- S  eV°T)  I<p(r+  sv0)|2  . 


(48) 


The  6  electric  constant  e(r>  is  related  to  the  index  of  refraction  by 

c<r>  =  nJ  . 

and  the  index  of  refraction  is  a  function  of  density  through  the  Lorentz-Lorenz  law 

n2  -  1 


n2  +  2 


=  Np  . 


(49) 


(50) 


In  Eq.  (50),  p  is  the  absolute  density  and  N  is  the  molecular  refractivity.  Since  p0  and 
p,  are  both  small  quantities,  Eq.  (50)  may  be  rewritten,  ignoring  quadratic  terms  in 
Pi  ,  as 


nJ  -  n0J  =  3Np,  . 

With  Relations  (51)  and  (48)  included,  Eq.  (1)  becomes 


2ik 


<P 


(51) 


-  ---(y--12-a---  e-“*(p  J  ds  (l- 8eV°T)  ^  |cp(rt 


sv0)|  =  0 


(52) 


which  is  the  fundamental  equation  to  be  solved  for  the  description  of  the  propagation  of  a 
laser  beam  in  the  presence  of  a  steady  wind.  The  boundary  conditions  involve  only  the 
specification  of  the  function  <p  on  the  plane  z  =  0.  Equation  (52)  is  clearly  nonlinear  and 
is  not  tractable  from  an  analytic  standpoint.  Therefore  it  is  solved  by  numerical 
computation. 


III.  COMPUTATIONAL  METHODS 
Transformation  to  Scaled  Variables 

Equation  (52)  is  scaled  for  purposes  of  numerical  computations.  Let  a  be  a  length 
characterizing  the  initial  beam  profile  in  a  direction  transverse  to  the  beam  axis.  New 
coordinates  x,  y,  and  S  are  defined  by 

x  =  x/a 

y  =  y/a  (53) 

£  =  z/kaJ  . 


NRL  REPORT  7293 


13 


A  dimensionless  wave  function  f  is  introduced  by  putting 


f(x.y.O.= 


/ca2 

y^*(x-y>z) 


(54) 


where  P  is  the  total  luminous  power  output  of  the  laser.  The  constant  of  the  motion, 
given  by  Eq.  (2),  can  be  stated  in  terms  of  f  by  the  requirement  that 


JJ  |f(x,y,£)|J  dx  dy  =  1  .  (55) 

C“ constant 

Combining  Eqs.  (53)  and  (54)  with  Eq.  (52),  the  equation  for  the  scaled  amplitude  f  is 

2i  —  +  -  /3e"^ka2“^  f  f  ds  (l  -  $eV°T)  |  f  |2  =  0  (56) 

where 

fi  =  3N(rc  V""  T7  '  (57) 

(When  a  slewing  beam  is  being  studied,  v0  in  Eqs.  (56)  and  (57)  is  to  be  replaced  by 
v0  +  fiz.) 

Because  Eq.  (56)  is  parabolic,  specification  of  the  amplitude  at  the  plane  £  =  0  de¬ 
termines  the  function  everywhere  else;  in  particular,  the  values  of  f  at  any  point  in  a 
given  plane  £  =  constant  can  be  determined.  The  technique  to  be  used  for  a  numerical 
computation  of  solutions  of  Eq.  (56)  will  be  the  replacement  of  Eq.  (56)  by  an  appropriate 
difference  equation  which  will  be  used  to  solve  for  the  solution  at  successive  planes  sep¬ 
arated  by  distances  limited  only  by  considerations  of  accuracy  and  stability.  The  method 
of  proceeding  from  plane  to  plane  is  referred  to  as  "marching."  (Because  of  the  char¬ 
acter  of  the  boundary  conditions  of  the  problem  and  of  the  difficulties  encountered  in  the 
numerical  solution  of  hyperbolic  equations  by  difference  techniques,  the  approximation  of 
neglecting  the  second  derivative  of  sp  with  respect  to  z,  used  to  derive  Eq.  (1),  is  essen¬ 
tial  to  a  numerical  solution  of  the  problem.  Marching  does  not  work  for  hyperbolic 
equations. ) 

Since  a  partial  differential  equation  does  not  uniquely  determine  a  corresponding 
difference  equation,  there  are  many  candidates  for  the  difference  equation  to  be  used  in 
the  numerical  computation.  Each  candidate  must  meet  the  requirements  of  convergence 
(to  insure  that  the  solution  of  the  difference  equation  converges  to  the  solution  of  the  dif¬ 
ferential  equation  for  any  £  in  the  limit  that  Ax  and  Ay  vanish)  and  of  stability  (which 
assures,  to  a  limited  degree,  that  the  numerical  solution  of  the  difference  equation  re¬ 
mains  reasonably  close  to  the  exact  solution  of  the  difference  equation.) 

Harmuth  (15)  has  shown  that  a  simple  forward  differencing  of  the  linearized  version 
of  Eq.  (56)  (obtained  by  replacing  the  integral  by  its  upper  bound)  is  unstable  for  all 
choices  of  step  size  A£  and  mesh  sizes  Ax,  Ay.  A  symmetric  differencing,  or  "two-point 
predictor,"  is  stable,  however,  for  a  nonempty  domain  of  the  parameters  A£,  Ax,  Ay. 

For  integration  in  £,  therefore,  the  derivative  in  £  is  replaced  by 
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II 

3C 


J.y.C 


f  (x,y.  5  +  AC)  -  f (x.y, £  -  A£) 
2A£ 


(58) 


while  the  algorithm  for  the  linearized  version  of  Eq.  (56), 


3f  32f  a2f 

2 i  — —  +  -  +  -  +  £Mf  =  0  , 

3C  3x2  3y2 


(59) 


is  taken  to  be 

..  [ftx.y.C  +  AC)  -  f(x.y.C-AC) 


f (x  +  Ax.y, C)  -  2f(x,y,C)  +  f ( x  -  Ax , y , C)"l 


(60) 

In  Eq.  (59),  m  is  taken  to  be  the  maximum  value  of  the  integral  in  Eq.  (56). 

The  stability  and  convergence  properties  of  the  linear  algorithm  given  in  Eq.  (60), 
which  are  necessary  (but  not  sufficient)  conditions  to  the  stability  and  convergence  of  the 
nonlinear  algorithm,  may  be  studied  (16)  by  examining  the  Fourier  components  of  all 
solutions  to  Eqs.  (59)  and  (60).  If  a  plane  wave  solution  to  Eq.  (59)  is  taken  as 

f  =  U(C)  (61) 

then 

U(C)  =  e-><p2  +  q2-/WK/2  (62) 

Assuming  a  solution  of  the  form  given  by  Eq.  (61)  for  the  difference  Eq.  (60),  the  u(c) 
now  has  two  solutions  (as  opposed  to  one  for  the  differential  equation): 

Uj(C)  =  ( iB  i  Vl  -  Bj)A(/2  (63) 


L 

f(x.y  +  Ay.Q  -  2f(x,y ,  Q  *  f(x,y-Ay,£) 

(Ay)2 

-AMffx.y.O  . 


where 


B  =  2AC 


sin2(pAx/2) 

(Ax)2 


s  in2(  qAy/2) 

(Ay)2  +  T ) 


It  is  necessary  to  show  that  there  exists  a  linear  combination  V(C)  of  the  two  solutions 
given  in  Eq.  f63)  that  converges  to  Eq.  (62)  in  the  limit  where  AC,  Ax,  and  Ay  vanish  in 
some  specified  order  and,  to  insure  stability,  that  V(C)  has  the  property  that  it  grows  as 


1  im  V  (  C)  ~  0 


(64) 
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If  (Ax)2  and  (Ay)2  are  required  to  vanish  in  proportion  to  Ac,  then  it  can  be  shown  that 


lim  U,(C)  =  c -HS'J-Wn-'i 
At-*o 


while  the  solution  U_(C)  fails  to  converge.  The  convergence  of  the  solution  of  the  differ¬ 
ence  equation  to  that  of  the  differential  equation  can  be  assured  if  there  exists  a  means 
of  suppressing  all  those  Fourier  components  of  the  type  U_(c).  The  unwanted  solution  in 
the  general  case  may  be  suppressed  by  a  "careful"  specification  of  the  function  f  at  the 
plane  C  =  0  and  at  the  plane  5  =  A 5.  By  "careful,"  one  means  that  the  specification  at 
5  =  A 5  must  be  at  least  as  accurate  (17)  as  a  Taylor's  expansion  about  5  =  0  to  the  same 
order  as  the  predictor  algorithm.  In  this  case,  the  algorithm  is  accurate  to  order  (AC)2, 
as  is  demonstrated  that 


f(C»AQ  -  f(g-AC) 
2A5 


1 

2A5 


1 

2A5 


Bf  1  ■>  "d2( 

f ( C)  +  A5  ^  ♦  A  (A?)2  —  ♦  0(A53) 


f(»  -  A5  ♦  i  (AO2  ~  ♦  6>*(AC3) 
°C  * 


=  0"(A53)  • 


(65) 


Thus  at  5  =  AC,  the  data  must  be  specified  as 


Bf  (AC) 

f  (AC)  -  f(0)  ♦  AC  (0)  t  — 

B2f 

5?<0) 

(66) 

where 

2i  (0)  =  -V2  f(x.y.O)  -  ftf  (x,y,0) 

J"  d.(|.s.^)  1  f  I1 
-»  ![0 

(67) 

and  B2f  (0)/BC2  must  be  similarly  specified. 

Since  the  unwanted  solution  has  been  suppressed  by  the  above  device,  the  test  for 
linear  stability  need  only  be  applied  to  the  wanted  solution,  and,  indeed,  may  be  applied 
to  its  Fourier  components  which  are  primarily  characterized  by  the  function  u,  (5).  The 
stability  requirement  on  ut(C)  given  by  Eq.  (64)  demands  of  B  that 


|  iB  +  Vl  -  B2i  <  1 


(68) 


Since  B  is  a  real  number,  if  B  >  1,  Eq,  (68)  can  never  hold.  A  marginal  stability  is 
achieved  only  if  B  <  1.  The  conditions  of  stability,  in  terms  of  AC,  Ax,  and  Ay,  divide 
into  three  categories: 


/3M  <  0  : 


0 


1  (±_  j_\  <  m  .  ac/3m  <  1 

2  \Ax2  A^2 )  4  ’  4  2 


l/3M|\  1 

4  /  2 


(69) 
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With  a  prescribed  choice  of  Ax  and  Ay  the  numerical  solution  of  the  difference  equation 
will  be  unstable  unless  AC  is  chosen  to  satisfy  the  three  conditions  in  Eq.  (69).  Since  the 
first  condition  is  the  most  stringent  of  all  three  restrictions,  choosing  AC  in  accord  with 
it  will  ensure  the  desired  stability. 

The  intrusion  of  PM  into  the  stability  conditions  given  in  Eq.  (69)  may  be  viewed  in  a 
different  light.  The  solution  of  Eq.  (59)  may  be  written  in  the  form 

£MC 

f  =  o  2  g 

where  g  is  a  solution  of  the  equation 

2i  +  V2g  =  0  . 

In  order  that  a  numerical  solution  of  the  equation  be  satisfactory,  the  oscillations,  in  C, 
of  the  exponential  factor  must  be  sampled  by  at  least  six  points  per  oscillation.  Thus, 
step  sizes  AC  must  be  chosen  so  that 


AC/3M  <  2 77/ 3  ~  2 


which  is  the  same  as  the  last  condition  of  Eq.  (69). 


All  the  above  considerations  apply,  of  course,  only  to  the  linear  approximation.  Not¬ 
withstanding  the  choice  of  M  as  the  maximum  of  the  integral,  the  predictor  was  found  to 
be  incapable  of  suppressing  a  nonlinear  instability  for  values  of  p  of  interest  even  though 
the  first  criterion  of  Eq.  (69)  was  well  satisfied.  In  order  to  damp  the  nonlinear  growth, 
a  simple  two-point  corrector  (18)  was  applied  in  conjunction  with  the  predictor.  Writing 
out  the  terms  of  order  AC3  explicitly,  the  predictor  and  corrector  are  taken  as 


f<p\C  +  AC)  =  f(t-AC)  ♦  2AC 


3f(0 

3C 


(AC)' 

3 


33f(C) 


and 


f(c)a  +  Ao  =  f(o  ♦  \  as 


9f<P>(C+AC) 

sc 


Sf(  0 

sc 


(AQ3  S3f(  0 
!2  ^ 


(70) 


When  AC  is  taken  small  enough,  the  third  derivatives  in  both  of  the  expressions  of  Eq. 
(70)  become  equal  in  value;  the  correct  value  of  f  at  the  plane  c f  AC  is  then  taken  to  be 
that  linear  combination  of  the  corrector  and  predictor  that  eliminates  the  terms  of  order 
AC3,  that  is, 


f  ( C  ^  AC)  -=  0.2f<P)(C  ♦  AC)  t  0.8f(c>(C  +  AC)  .  (71) 

In  Eq.  (70),  the  first  derivatives  are  evaluated  by  the  use  of  the  algorithm  that  relates 
these  derivatives  to  the  transverse  derivatives  and  the  heating  integral  at  the  plane  C . 
just  as  Eq.  (67)  does  so  at  the  plane  which  is  the  face  of  the  laser.  This  form  of  the  algo¬ 
rithm  has  been  found  to  be  sufficient  to  suppress  the  nonlinear  growth  and  has  been  used 
in  this  work. 

It  is  convenient  (19)  to  make  a  change  of  variables  in  order  to  map  the  whole  trans¬ 
verse  plane  into  a  finite  region.  The  following  transformations  were  chosen: 


i 
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x  =  In  [#'(4-#)] 

y  =  In  [t)/(4  -  7))]  . 

(72) 

The  transverse  Laplacian,  in  terms  of  the  new  variables  £  and  t),  now  becomes 

£2(4-f2)  3Jf  .  (4-2f)(4-£)  3f 
16  +  16  W 

V2(4- v1)  3Jf  (4  -  2tj)C4  -  tj) 

T  16  +  16  ^  ‘ 

In  terms  of  these  new  variables,  the  conserved  quantity  (see  Eq.  (54))  takes  the  form 


d-r) 


f  *  *  f* 


(4 -f)tf(4 -■»)!? 


1  . 


(74) 


The  left-hand  side  of  Eq.  (74)  is  monitored  throughout  every  run  as  a  check  on  the  sta¬ 
bility  and  convergence  of  the  algorithm. 

Laser  beams  whose  transverse  dimensions  vary  considerably  with  z  present  prob¬ 
lems  of  sampling  for  numerical  schemes  with  fixed  transverse  mesh  size.  In  Appendix 
A,  a  formulation  of  the  propagation  problem  is  presented  in  terms  of  coordinates  that 
will  automatically  take  into  account  the  change  in  the  beam  size.  The  stability  and  con¬ 
vergence  criteria  developed  in  this  section  may  be  applied  to  the  equations  of  Appendix  A 
with  little  modification. 


IV.  RESULTS  OF  NUMERICAL  COMPUTATIONS 
Accuracy  Criteria 

The  computer  code  described  in  Section  III  has  been  run  successfully  for  a  large 
number  of  initial  parameters.  Because  exact  analytical  solutions  of  the  full  nonlinear 
problems  were  lacking,  criteria  were  needed  for  determining  whether  a  given  run  led  to 
correct  results  or  not.  Since  comparisons  with  the  exact  solutions  could  be  the  only 
valid  check  on  the  numerical  solutions,  the  computer  code  was  applied  to  problems  to 
which  exact  results  were  known.  These  solutions  are  discussed  in  Appendix  B.  Such 
checks  provided  a  testing  ground  for  certain  parts  of  the  numerical  procedures  used  and 
served  as  a  test  bed  for  studying  the  effects  of  parameter  such  changes  as  step  sizes, 
mesh  sizes,  focused  coordinate  systems,  etc.  These  comparisons  provided  the  necessary 
minimum  level  of  confidence  in  the  chosen  algorithm  needed  to  proceed  with  the  nonlinear 
problems. 

This  section  is  devoted  to  a  description  of  the  criteria  that  were  and  are  used  to 
establish  run  quality,  the  characteristics  of  reliable  and  unreliable  runs,  and  a  delinea¬ 
tion  of  the  potential  sources  of  error.  These  are  illustrated  by  the  presentation  of  some 
of  the  results  obtained  with  the  computer  code. 
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Run  Quality  Criteria 

Run  quality  is  a  function  of  the  downrange  distance  z;  a  given  run  may  be  good  up  to 
one  value  of  z ,  then  poor  thereafter.  Clearly  it  is  desirable  to  have  runs  of  high  quality 
for  as  large  a  value  of  range  as  possible. 

It  has  been  shown  that  the  quantity 


is  a  constant  of  the  motion.  The  ratio 


SE(z)  /E  = 


E<«)  -  E(0) 
E(0) 


provides  a  quantitative  measure  of  numerical  accuracy  in  numerical  solutions.  This 
quantity  was  determined,  for  the  runs  to  be  described  later  in  this  section,  at  selected 
intervals  downrange,  but  in  this  report,  its  value  will  be  given  only  at  the  last  oife  or  two 
values  of  ^  before  run  termination.  Since  se/e  generally  (but  not  always)  increases 
with  z,  the  final  numbers  are  indicative  of  the  accuracy  upbeam. 

A  second  qualitative  measure  of  the  validity  of  a  calculation  is  the  frequency  with 
which  the  oscillations  in  the  real  and  imaginary  parts  of  the  amplitude  are  sampled  by 
the  mesh.  To  each  run,  usually  at  the  final  one  or  two  z  values  reported,  a  Run  Quality 
Factor  (RQF)  is  assigned  according  to  the  following  conditions: 


Sampling  Characteristics 

Sampling  nowhere  exceeds  6  to  8  points  of 
the  mesh  per  oscillation. 

Small  amplitude  oscillations  sampled 
poorly;  large  amplitude  oscillations  sam¬ 
pled  at  least  6  to  8  times. 

Large  and  small  oscillations  are  sampled 
at  least  6  to  8  times. 


Run  Quality  Factor 
POOR 


BORDERLINE 


GOOD 


Large  and  small  oscillations  are  sampled 
better  than  8  times. 


EXCELLENT 


Figures  2  and  3  show  a  plot  of  the  real  part  of  the  amplitude  at  the  mesh  points  taken 
from  a  typical  set  of  runs  demonstrating  the  sampling  for  different  RQF's.  Run  quality 
factors  are  helpful  in  discussing  run  characteristics,  but  they  must  be  used  with  caution 
since  they  can  be  misleading.  An  RQF  of  POOR  may  in  fact  be  describing  the  oscillations 
quite  faithfully;  but  usually  poor  runs,  when  continued  downrange,  lead  rapidly  to  very 
large  values  of  £*'/e.  A  run  with  an  RQF  of  EXCELLENT  may  be  suffering  from  aliasing 
and,  in  fact,  be  quite  poor;  cases  such  as  these  may  often  be  detected  by  looking  at  RQF 
as  a  function  of  i . 


Further,  RQF's  are  often  difficult  to  assign;  they  are  a  function  of  mesh  size  as 
well,  and  they  suP”*-  a  degree  from  subjective  interpretations.  A  cruder  pair  of 
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BORDERLINE 


ir 

POOR 


I-  u  m> 


6106;  P  =M06  watts;  a  =  25  cm;  f  *  1.0  km;  L  =1.2  km; 

•O.  =  .05  rad/sec;  p^O)  =  12  torr;  v0  =  1000  cm/sec. 

Fig.  2  -  Sampling  quality  of  the  real  part  of  the  amplitude  along  the  line  of  sym¬ 
metry  graphically  displayed  for  Run  Quality  Factors  BORDERLINE  and  POOR. 


categories  such  as  acceptable  runs  and  unacceptable  runs  may  in  some  instances  be  sat¬ 
isfactory.  Experience  shows  that  POOR  runs  are  unacceptable ,  while  BORDERLINE, 
GOOD,  and  EXCELLENT  runs  are  acceptable . 

A  third  and  very  important  criterion  of  run  quality  is  provided  by  the  appearance  of 
isoirradiance  contours  in  a  beam  cross  section  as  a  function  of  z .  For  a  run  of  high 
quality,  these  contours  evolve  smoothly  from  concentric  circles  (for  an  initially  rota- 
tionally  symmetric  beam)  to  crescent-shaped  curves.  Runs  of  poor  reliability  generally 
betray  themselves  by  the  development  of  straight  lines  about  which  the  crescent  patterns 
are  forced  to  evolve.  These  lines,  which  are  not  isoirradiance  contours,  and  hence  are 
not  plotted,  run  parallel  to  the  mesh  lines  or  at  45°  to  the  mesh  lines.  There  are  two 
principle  reasons  for  such  development.  First,  in  the  case  of  a  small  f  ratio,  the  initial 
amplitude  will  contain  many  oscillations  in  the  aperture  plane.  A  fine  mesh  is  needed  to 
sample  these  sufficiently  frequently  for  a  good  description.  Since  a  rectangular  mesh  is 
used,  sampling  along  the  45°  lines  is  less  frequent  than  along  the  mesh  lines  themselves. 


20 


AITKEN,  HAYES,  AND  ULRICH 


ir 

yr.'TT.T.BIT 


z 


0 


Fig.  3  -  Sampling  quality  of  the  real  part  of  the  amplitude  graphically 
displayed  for  Run  Quality  Factors  EXCELLENT  and  GOOD. 


Mesh  sizes  that  provide  adequate  sampling  at  the  laser  aperture  may  prove  inadequate  at 
some  distances  downbeam,  and  errors  are  introduced  into  the  numerical  solution. 

In  the  vacuum  case,  the  errors  possess  a  90°  rotational  symmetry.  When  wind  is 
present,  the  symmetry  in  the  errors  deteriorates  and  is  masked  to  some  degree  by  the 
aberrations  due  to  heating.  These  effects  are  illustrated  by  the  contour  plots  of  Figs. 

4  and  5.  Figure  4  shows  two  vacuum  runs  for  a  Gaussian  beam  focused  at  1  km;  the 
mesh  number  was  31  x  61  in  both  cases,  but  the  sampling  of  the  oscillations  was  altered 
by  using  two  different  coordinate  systems.*  The  contours  in  the  top  row  were  calculated 
with  i  =  1.5  km,  while  those  of  the  bottom  row  were  calculated  with  t  =  1.12  km.  In  the 
first  case,  the  beam  is  sampled  well  enough  by  the  mesh  at  the  laser  aperture,  but  as  the 
beam  converges  toward  the  focus,  the  sampling  becomes  poorer.  In  the  second  case, 
improvement  is  achieved  for  two  reasons: 


*See  Appendix  A  for  a  discussion  of  focused  coordinate  systems  and  the  significance  of  the  param¬ 
eter  2  mentioned  in  the  above  paragraph. 


beams  focused 
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1.  The  beam  is  being  specified  on  spherical  surfaces,  on  which  phase  oscillations 
are  small  since  the  surfaces  more  closely  approximate  surfaces  of  constant  phase  than 
planes  do. 

2.  The  mesh  size  shrinks  in  about  the  same  proportion  that  the  beam  does,  so  that 
the  beam  is  sampled  with  the  same  thoroughness  downrange  as  it  is  at  the  aperture. 

Figure  5  shows  a  run  with  atmospheric  heating  present,  but  with  an  inappropriate 
choice  of  focused  coordinate  system.  The  outermost  contours  betray  the  poor  sampling 
on  the  left  (the  wind  moves  from  left  to  right),  but  the  aberrations  of  the  beam  due  to 
atmospheric  heating  mask  this  effect  on  the  downwind  side.  Secondly,  contours  may  show 
an  evolution  about  lines  parallel  to  the  mesh  axes;  such  development  arises,  in  part,  for 
reasons  of  physics.  Any  numerical  scheme  for  solving  Eq.  (56)  must  rely  on  a  mesh 
which  samples  only  a  part  of  space;  a  beam  initially  confined  to  that  part  of  space  cov¬ 
ered  by  the  mesh  will  ultimately  grow  in  size  unless  self-trapping  should  occur.  This 
growth  occurs  both  for  vacuum  beams  and  for  beams  traversing  a  medium.  Unless  the 
mesh  is  chosen  so  as  to  be  able  to  accommodate  this  increase  in  size,  a  point  will  be 
reached  in  the  calculation  where  the  beam  amplitude  in  the  outer  regions  of  the  mesh  is 
large,  and  poor  sampling  or  boundary  conditions  will  begin  to  manifest  themselves  in  the 
solution.  The  isoirradiance  contours  become  very  complicated  and  show  definite  square 
edges,  reflecting  the  presence  of  the  edges  of  the  mesh;  for  this  reason,  we  describe 
such  a  situation  by  the  phrase  "the  beam  has  hit  the  edges  of  the  mesh."  Run  quality  at 
such  downrange  distances  deteriorates  completely,  reflecting  poor  sampling,  and  insta¬ 
bilities  arise  leading  to  large  values  for  se/e.  Figure  7  shows  the  isoirradiance  con¬ 
tours  for  such  a  case,  which  will  be  discussed  below  in  greater  detail.  In  some  instances, 
the  options  made  available  by  the  parameter  8  in  the  focused  coordinate  systems  enable 
the  deterioration  of  the  run  to  be  postponed  to  larger  values  of  z.  However,  the  condi¬ 
tions  imposed  upon  the  value  of  8  to  improve  sampling  at  the  aperture  for  beams  of  small 
f  ratio  run  counter  to  those  arising  from  the  need  to  avoid  having  the  beam  hit  the  edges 
of  the  mesh.  In  cases  where  neither  condition  may  be  relaxed,  the  computer  code  cannot 
be  relied  upon  to  describe  the  correct  physical  situation.  Modifications  of  the  code  that 
will  allow  8  to  vary  with  z  in  an  appropriate  manner  appear  to  be  a  way  of  handling 
situations  such  as  these,  but  such  changes  have  not  yet  been  incorporated  into  the 
program. 


Numerical  Results 

Results  of  particular  numerical  solutions  of  the  problem  of  laser  beam  propagation 
through  a  nonturbulent  atmosphere  are  now  presented.  The  numbers  attached  to  a  given 
computation  serve  to  designate  it;  the  particular  characteristics  of  the  beam  at  the  aper¬ 
ture  are  specified,  followed  by  the  run  quality  characteristics  (5E/fe  and  RQF)  at  one  or 
more  distances  downrange.  The  figures  show  the  isoirradiance  contours  at  selected 
downrange  distances.  (In  these  figures,  the  scale  of  each  of  the  plots  is  set  by  the  size 
of  the  first  of  the  plots;  the  contours  represent  90%  of  peak  intensity,  80%  of  peak  inten¬ 
sity,  etc.,  from  the  inside  out.  The  actual  value  of  the  peak  intensity  and  the  deflection 
of  the  peak  intensity  point  from  the  beam  axis  are  data  available  from  the  computations 
in  general,  but  in  some  of  these  early  results  shown  here,  they  were  not.)  A  paragraph 
is  devoted  to  a  discussion  of  the  results  in  each  case.  It  is  emphasized  here,  and  in  the 
examination  of  individual  runs,  that  not  all  of  the  solutions  presented  here  are  to  be 
regarded  as  acceptable.  Indeed,  unacceptable  numerical  results  are  presented  to  illus¬ 
trate  the  characteristics  of  erroneous  solutions  that  may  arise  from  an  uncritical  use  of 
the  numerical  procedures  described  in  this  report. 
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FOCUSSED  FhOPAGAl ION  IN  CONSTANT  INDEX  MEDIUM 


Fig.  6  -  Comparison  of  computer  results  with  analytical  results  obtained  from  the  theory  of 
propagation  of  a  focused  laser  beam  in  vacuum;  a  =  10  cm. 


The  graphs  in  Fig.  6  show  light  rays  emanating  from  a  point  at  the  laser  aperture 
where  the  power  is  30%  of  the  central  intensity.  The  runs  here  are  in  vacuum  to  check 
the  linear  aspects  of  the  code.  The  solid  curves  are  analytical  results,  the  points  rep¬ 
resent  the  computer  results,  and  the  arrows  point  to  the  waist  of  the  beam.  The  agree¬ 
ment  between  theory  and  numerical  computation  is  excellent.  The  runs  in  Fig.  6  were 
made  with  8  =  1.5  f  in  all  cases.  Because  the  beam  diverges  to  the  right  of  the  waist 
points,  while  the  focused  coordinate  system  continued  to  converge,  the  beams  ultimately 
hit  the  edges  of  the  mesh.  The  quantities  SE/E,  at  the  focal  distances  for  the  cases 
shown,  were  6xlCT4,  7  xlO'4,  8xl0"4,  and  2xl0'3  for  the  cases  f  =  1,  2,  3,  and  4  km, 
respectively,  while  for  f  =  5  and  f  =  6  km,  SE/E  was  2  xlO'3  just  short  of  the  focal  point 
in  each  case. 


No-Diffraction  Runs 

No  figure  is  provided  here.  The  rationale  for  these  runs  is  provided  in  Appendix  B. 
These  runs  are  made  to  check  the  nonlinear  aspects  of  the  code.  This  is  done  by  drop¬ 
ping  the  Laplacian  term  in  the  wave  equation,  which  is  equivalent  to  saying  that  refrac¬ 
tion  dominates  the  character  of  the  beam  propagation.  Such  runs  might  be  termed  high- 
power  runs,  or  no-diffraction  runs,  or  the  geometrical  optics  limit.  The  resultant 
equation,  while  still  nonlinear,  may  be  solved  exactly  (see  Appendix  B)  to  give  the  result 
that  the  intensity  profile  does  not  change  downbeam.  This  result  was  confirmed  by 
computation  with  the  code  by  circumventing  the  subroutines  that  compute  the  Laplacian. 


If  raw 
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Run  No.  3346,  shown  in  Fig.  7,  has  the  same  parameters  as  those  of  the  next  two 
runs,  illustrated  in  Figs.  8  and  9.  This  sequence  of  runs,  together  with  the  vac  .urn  run 
for  2  km,  shown  in  Fig.  6,  illustrates  the  potential  sources  of  error  that  may  be  encoun¬ 
tered  in  these  numerical  computations  as  discussed  earlier  in  this  section.  With  no  prior 
knowledge  of  the  extent  to  which  the  beam  would  bloom,  or  diverge,  because  of  the  heat¬ 
ing,  a  first  guess  was  made  that  the  choice  of  focused  coordinates  used  for  the  vacuum 
case  would  work  here.  The  first  indication  that  this  assumption  was  erroneous  was  pro¬ 
vided  by  the  isoirradiance  contours  at  z  =  1.0098  km;  the  maximum  linear  distance  of 
the  0.3  isophote  from  the  horizontal  center  line  is  the  same  as  it  was  at  the  laser  aper¬ 
ture  while,  for  a  vacuum  run  as  seen  from  Fig.  6,  it  is  only  one-half  as  large.  Thus,  the 
heating  effect  is  preventing  the  focusing  of  the  beam,  as  it  would  otherwise  do  in  vacuum. 
The  coordinate  system,  on  the  other  hand,  is  still  converging  so  the  beam  will  approach 
the  edges  of  the  mesh.  At  z  =  1.2108  km,  discernible  flattening  of  the  contours  at  right 
angles  to  the  wind  direction  is  perceived  and  become  more  pronounced  at  z  =  1.4117  km. 
The  structure  at  this  point  cannot  be  believed,  even  though  SE/E  =  7  xlO'4.  (For  this 
run,  the  RQF  was  not  available.  A  reasonable  guess  can  be  made  from  the  rapid  varia¬ 
tions  of  intensity  across  the  beam  that  RQF  here  is  POOR. )  Succeeding  isophote  plots 
demonstrate  the  total  deterioration  of  the  run.  At  z  =  1.6127  km,  SE/E  ~  10 208,  while 
at  z  =  1.8337  km,  SE/E  =  5.91. 

It  is  important  to  note  that  the  patterns  at  z  =  1.6127  km,  which  so  clearly  show  the 
unacceptable  run  quality,  could  easily  have  been  missed  had  the  computer  been  com¬ 
manded  to  plot  at  slightly  different  z  values.  If  fewer  isoirradiance  contours  had  been 
used,  the  plots  at  z  =  1.2608  km  and  1.4117  km  would  have  shown  less  complexity  and 
could  then  have  been  readily  interpreted  as  complicated  diffraction  patterns.  However, 
such  interpretations  were  avoided  by  use  of  the  SE/E  and  RQF's. 

The  initial  data  shown  in  Fig.  8  are  the  same  as  the  previous  run;  '\owever,  £  =  a>t 
i.e.,  ordinary  Cartesian  coordinates  were  used,  while  (Nx,Ny)  =  (31,61). 

Because  the  previous  rim  was  so  clearly  in  error,  it  became  necessary  to  do  the 
computation  again.  Since  the  beam  did  not  appear  to  be  focusing,  Cartesian  coordinates 
were  used.  At  z  =  2  km,  se/e  =  7  xlO"5,  but  no  RQF's  were  available  in  this  early  run. 
However  the  contours  at  z  =  2.0  km  show  considerable  flattening  on  the  wind  side  of  the 
beam.  Beam  deflection  into  the  wind  was  enough  to  cause  concern  about  poor  sampling, 
so  it  was  therefore  decided  that  the  run  was  UNACCEPTABLE. 

Run  No.  3353  of  Fig.  9  is  a  repeat  of  the  above  two  runs,  but  with  £  =  -2  km.  This 
choice  of  £  was  taken  because  the  coordinate  system  would  then  be  slightly  diverging, 
giving  the  beam  more  space  in  which  to  diverge.  The  differences  between  the  isoir¬ 
radiance  contours  at  z  =  2  km  in  this  case  and  in  the  previous  run  are  quite  pronounced. 
The  squaring  of  the  pattern  as  seen  in  Fig.  8  is  gone,  and  the  contours  on  the  wind  side 
■*re  less  compressed.  The  value  SE/E  =  2  x  10‘ 5  here,  and  the  RQF  was  GOOD.  It  ap¬ 
pears  that  the  run  could  be  further  improved  by  choosing  £  slightly  closer  to  zero;  how¬ 
ever,  this  would  sacrifice  accuracy  at  the  center  of  the  beam  for  the  same  total  number 
of  mesh  points.  With  the  NRL  computer,  increasing  the  number  of  mesh  points  is  not 
possible. 

The  contours  show  that  the  beam  reaches  a  minimum  size  at  approximately  1.2  km 
and  diverges  for  larger  values  of  z.  Since  the  computer  printout  also  included  phase 
angle  at  each  mesh  point,  the  assessment  that  the  beam  is  no  longer  focusing  beyond  1.2 
km  could  be  verified  by  the  behavior  of  the  phase  angle  as  a  function  of  radial  distance 
from  the  beam  center.  Phase  angle  should  increase  with  radial  distance  for  a  diverging 
beam,  and  this  was  indeed  the  case  here. 
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Fig.  9  -  Run  No.  3353. 


The  deflection  of  the  peak  intensity  along  the  line  of  symmetry  of  the  contour  plot  is 
plotted  against  range  z  in  Fig.  10.  Deflection  at  small  z  grows  quadratically  with  z ,  as 
geometrical  optics  predicts.  At  1.8  km,  the  deflection  curve  show  a  rather  severe  change 
in  slope.  This  appears  to  be  due  to  a  rapidly  changing  shape  in  the  profile  and  to  the  fact 
that  the  beam  is  sampled  a  finite  number  of  times. 

Beam  deflection,  measured  in  the  above  fashion,  is  a  crude  measure  of  the  location 
of  the  beam.  It  is  quite  conceivable  that,  with  the  growth  of  diffraction  peaks  along  the 
line  of  symmetry,  the  deflection  itself  can  show  an  abrupt  discontinuity  while  the  beam 
configuration  is  changing  in  a  smooth  and  continuous  fashion.  Therefore  the  position  of 
the  peak  intensity  point  along  the  central  line  indicates  the  beam  displacement  in  a 
qualitative  sense  only. 
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RANGE  (KM) 

Fig.  10  -  Beam  deflection  vs  downrange  distance  for  Run  No.  3353. 


The  value  of  the  peak  intensity  along  the  line  of  symmetry  in  a  cross  section  of  the 
beam  is  plotted  vs  range  z,  in  Fig.  11.  The  peak  intensity  is  seen  to  reach  a  maximum 
at  1.2  km.  The  decrease  thereafter  is  related  directly  to  the  divergence  of  the  beam. 

The  isoirradiance  contours  of  Run  #1155  are  shown  in  Fig.  12.  At  z  =  1  km, 

6E/E  =  1  x  10"  3  and  the  RQF  is  EXCELLENT.  The  sharp  corners  of  the  interior  contours 
can  be  eliminated  by  using  a  finer  mesh.  Slewing  and  kinetic  cooling  are  present.  At  a 
range  of  1  km,  v0r  =  30  cm,  while  the  beam  size  is  reduced  to  about  6  to  10  cm  in  diam¬ 
eter.  Thus,  even  with  a  large  amount  of  water  vapor  present,  the  air  is  swept  out  of  the 
beam  before  the  excited  nitrogen  molecules  decay  to  add  to  the  heating,  reducing  the 
overall  heating  effect  by  .6.  The  high  slewing  rate  further  diminishes  the  heating.  Al¬ 
though  the  calculation  was  not  carried  further,  it  ;s  clear  that  diffraction  will  cause  the 
beam  to  spread  at  larger  ranges. 

Beam  deflections  were  less  than  1  cm  at  all  ranges  up  to  1  km  in  this  run. 

Plotted  in  Fig.  13  are  the  isoirradiance  contours  at  the  focal  point  for  runs  which 
have  identical  parameters,  except  water  vapor  pressure,  to  demonstrate  the  effect  on  a 
focused  beam  of  the  kinetic  cooling  phenomenon.  All  figures  are  drawn  to  the  same  scale. 
An  absolute  scale  is  shown  below  the  left-most  figure.  The  peak  intensities,  from  left  to 
right,  are  3140  watts/cm2,  490  watts/cm2,  and  640  watts/cm2.  The  peak  intensity  value 
for  the  central  figure  is  lower  than  that  of  the  right  one  because  of  the  development  of 
two  strong  maxima  in  the  beam.  In  contrast,  the  1/e  radius  of  a  diffraction-limited  beam 
(shown  as  a  dashed  circle  centered  on  the  beam  axis,  indicated  by  a  small  cross)  is  1.65 
cm,  with  a  peak  intensity  of  11,250  watts/cm2. 


MAXIMUM  INTENSITY  (W/CM«) 


.  12  -  Run  No.  1155. 
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Fig.  13  -  Comparison  of  the  effects  of  water  vapor  pressure  P(H20)  on  a  focused  beam. 


\ 


The  parameters  that  characterize  the  above  runs,  apart  from  water  vapor,  are 
P  =  10s  watts;  a  =  10  cm;  v0  =  200  cm/sec;  fi  =  0.00  rad/sec;  f  =  1  km;  and  (n  ,  N  )  =  I 

(31,61). 

Run  #3145,  whose  contours  are  shown  in  Fig.  14,  was  chosen  with  parameters  to 
match  the  run  reported  by  Wallace  and  Camac  (11),  who  studied  laser  propagation  using  < 

geometrical  optics.  They  presented  their  contour  plot  at  z  =  .6  km.  The  two  figures, 
with  the  different  normalizations  accounted  for,  are  in  detailed  agreement.  Beyond  this 
point,  at  798  km,  the  beam,  which  was  initially  collimated,  is  undergoing  self-focusing. 

(At  z  =  798  m,  'E/E  =  3x  10“4  and  RQF  is  EXCELLENT.)  The  next  figure  can  only  be 
regarded  as  a  qualitative  representation  of  the  phenomenon  because  the  beam  is  collaps¬ 
ing  to  too  small  a  configuration  to  be  handled  by  the  mesh  size  used  in  the  computation. 
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Appendix  A 

FOCUSED  COORDINATES 


In  numerical  computation  of  solutions  to  Eq.  (56),  sufficiently  small  mesh  sizes  to 
describe  the  oscillations  of  the  amplitude  in  the  transferse  directions  are  imperative. 

If,  for  focused  beams,  a  fixed  transverse  mesh  is  used,  there  may  be  insufficient  sam¬ 
pling  in  the  vicinity  of  the  focal  region,  especially  for  short  focal  lengths  where  the  spot 
size,  in  vacuum  at  least,  is  certainly  small.  Variable  mesh  sizes  may  be  achieved  by 
using  angular  coordinates  on  concentric  spherical  surfaces  as  follows: 

As  was  done  in  Section  III,  construct  a  right-handed  coordinate  system  whose  origin 
is  at  the  center  of  the  laser  face,  whose  plane  lies  in  the  xy  plane,  and  whose  ?.  axis  is 
the  direction  of  propagation.  Downbeam  a  distance  2  p’ece  the  origin  of  a  second  right- 
handed  coordinate  system,  with  its  z'  axis  pointing  back  towards  the  laser  (see  Fig.  Al). 


Fig.  Al  -  Relationship  between  three  coordinate 
systems  used  in  the  discussion  of  focused  coordi¬ 
nates. 


With  respect  to  this  system  construct  the  customary  spherical  polar  coordinates;  the 
coordinates  of  a  point  P  in  space  will  thus  be  characterized  by  (x,y,z)  or  by  (r',  e',  q>'). 
The  Helmholtz  equation  in  terms  of  the  latter  variables  is 


1  3  ,  ,.2  30  1  3  .  30  1  320  .  .  2  ,  n 

(r')J  *r*  3"  (r')2  sine  35'  30'  ( r  <  )2  s  in2<r  3((p«  )J  *  (A1) 


For  a  wave  propagating  along  the  positive  z  direction  we  put 
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where 


k  -  kej/2  +  ia/2  , 

into  Eq.  (Al)  and  derive  an  equation  for  <P  which  will  involve  32<t>/3(  r')2  and  k  3i>/3r'. 
Dropping  the  second  space  derivatives  as  was  done  in  the  derivation  of  Eq.  (1)  and  carry 
ing  only  lowest  order  terms  in  the  absorption  coefficient,  the  equation  for  <i>  becomes 


-2ikc 


1/2 

0 


1  3  .  3<t>  1  s2<k 

- sin#  —  4-  -  - - 

sin#'  3#'  30'  sin2#' 


+  k2(e<r>  -  e$r>)  <i>  =  0  . 


(A3) 


At  points  between  the  focal  point  f  of  the  beam  and  the  laser  face,  the  amplitude 
will  usually  be  small  for  all  but  small  values  of  #' .  Carrying  only  the  lowest  orders  of 
#',  Eq.  (A3)  is  closely  approximated  by 


-2ike 


1/2 

0 


3<t>  1 

‘i  J_  I 

V  Ji.) 

li  1  ^  1 

3r'  (rJ)2 

#'  30'  \ 

,  »■) 

'  (O')2  3(V')2 

+  k2(e<r>  -  e<r>)  <t  =  0  . 
A  coordinate  transformation 


converts  Eq.  (A4)  into 


x  -  O'  cos<p'  ,  y  =  #'  sin<p' 


(A4) 


(A5) 


-2ike 


1/  2 
0 


3r  ‘ 


i  f  32<t>  a2<i>\ 

(r')2\3x2  +  3p2/ 


k2(c<r>-  4r>)4>  =  0 


(A6) 


which  is  quite  similar  in  form  to  Eq.  (1).  The  quantities  x,  and  y  are  angular  coordi¬ 
nates  on  spheres  of  radius  r '  in  the  transverse  directions,  as  seen  from  the  origin  of 
the  primed  coordinate  system. 


If  a  is  a  characteristic  transverse  linear  dimension  (say  beam  radius)  characteriz¬ 
ing  the  initial  laser  beam  profile,  then  #0  =  -i/£  is  a  corresponding  characteristic 
angular  dimension.  Scaling  the  variable  x  and  y  with  this  quantity,  a  new  pair  of  varia¬ 
bles  x  and  y  are  introduced: 


X 


-  t  ~ 

y  =  —  y  • 

J  a  3 


(A7) 


For  those  points  in  space  where  the  amplitude  is  significantly  different  from  zero,  the 
variable  r'  may  be  related  to  the  propagation  distance  z,  to  lowest  order  in  ,  by 

r'  =  £  -  z 

and  z  is  related  to  {  by  Eq.  (53).  In  terms  of  the  variables  x,  y,  and  c»  Eq*  (A6) 
becomes 
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2i  *r 


while  the  intensity  I  now  takes  the  form 


Vi2<t  +  kV(c<0_£(r)>  ^  .  o 


(A8) 


I  =  0->tn2£l!, 

,2 


77a 


(1-kaVf)2 

The  index  of  refraction  is  related  to  the  amplitude  $  by 


(A9) 


-f-7) 


1*1 


(A10) 


The  variables  x,  y,  and  5  may  be  related  to  the  original  variables  x,  y,  and  z  by 


X=(1-|)ay,  y=(l-|)ax.  z=ka2C. 

In  the  limit  that  £ — ►<»,  all  previous  results  are  recovered. 


(All) 


The  solution  to  the  boundary  value  problem  posed  by  the  differential  equation  [Eq. 
(A9)]  requires  the  specification  of  the  function  on  the  spherical  surface  that  goes  through 
the  point  x  =  y  =  z=0.  On  such  a  surface,  the  amplitude  of  a  Gaussian  beam  focused  at 
a  distance  f  is  given  by 


1 


ik£ 


irr 


(M2) 


To  describe  the  oscillations  in  the  real  and  imaginary  parts  of  the  amplitude  that 
are  due  to  the  last  factor  in  Eq.  (A12),  a  transverse  mesh  size  of  appropriately  small 
dimensions  must  be  chosen.  The  larger  the  beam  dimensions,  the  smaller  this  mesh  size 
must  be;  the  mesh  size  will  also  be  correlated  with  the  f  ratio  of  the  beam.  For  values 
of  £  near  f,  such  oscillations  can  be  reduced  considerably;  for  £  =  f ,  they  vanish.  If  one 
is  not  inteiested  in  the  behaviour  of  the  beam  as  far  as  the  focal  point,  this  latter  choice 
of  £  is  the  appropriate  one.  Larger  values  of  £  are  required  to  study  the  beam  in  the 
vicinity  of  the  focal  length,  since  the  denominators  of  the  second  and  third  terms  (see  Eq. 
(A10))  in  the  wave  equation  (Eq.  (A8))  vanish  at  z  =  £,  and  the  approximations  on  which 
the  equation  was  based  are  no  longer  valid. 

The  numerical  value  of  the  quantity  £  is  determined  by  the  problem  itself.  It  is  con¬ 
venient,  in  practice,  to  attempt  a  numerical  solution  with  fixed  mesh  size.  If  the  results 
show  a  strong  focusing  of  the  beam,  the  numerical  results  can  be  made  more  reliable  by 
redoing  the  computation  with  £  chosen  slightly  larger  than  the  focal  length.  If  the  heat¬ 
ing  causes  severe  blooming  in  spite  of  the  initial  focusing,  it  may  prove  necessary  to  use 
a  negative  value  of  £  so  that  the  coordinates  are  defocusing. 

An  alternative  derivation  of  Eq.  (A9),  which  is  much  simpler  than  the  above  but  has 
less  immediate  geometrical  significance,  is  obtained  by  putting 


38 


AITKEN,  HAYES,  AND  ULRICH 


so  that 


g.ving  the  equation 


^  l  .  ,  %  l  _ 

x  -  - - -  x  and  y  =  - - - y 

1  -  ka2£/f  1  -  ka2^f 


f(x,y,L  =  f(x(x.O.  y(y, 0.  0 


3f  ..  1  3f  _3f\  1  ~2- 

+  l-ka2£/f  \  &  f  y  2y)  +  (l-kaVi)J  1 


+  k  2a2(  £  -  £0)  f  =  0  . 

The  first-derivative  terms  are  eliminated  by  defining  a  new  amplitude  <p  by 


.  ka2  /  ka2C\  2 

r -  cxpr  ^  \  r)(x 


1  -  ka2C/£ 


<J>(x,  y.O  • 


<f  is  then  the  amplitude  in  Eqs.  (A9)-(A11).  The  equivalence  of  the  two  derivations  is 
clear. 


Appendix  B 

ANALYTICAL  RESULTS 


The  computer  code  can  be  checked,  in  part,  by  comparing  computer  solutions  of 
known  problems  with  analytical  results.  The  analytical  results  presented  here  provide 
a  check  on  two  different  aspects  of  the  code. 


Vacuum  Propagation 

The  vacuum  propagation  of  a  laser  beam  is  described  by  Eq.  (1)  with  /3  set  equal  to 
zero.  The  solution  of  the  resultant  equation  is 


<p!p,  z)  = 


=  -  k  t JJ,lx'(,y'  ) 


(Bl) 


The  integration  is  carried  over  the  plane  z  =  0,  and  <p(p,0)  is  the  amplitude  specified  at 
the  face  of  the  laser.  For  a  Gaussian  beam  whose  power  is  normalized  to  unity  and 
whose  intensity  is  reduced  from  its  central  value  by  a  factor  of  e  at  a  radial  distance  a, 
and  which  is  focused  at  a  distance  f  down  the  beam  axis, 


V(P.O)  — rrexp 

AT2 


L-  • 

2  ‘  2f 


and 


,  -i/  /71a2 

"(p'° 1  rrmTn oxp 


a 


2[C2  -  f/f)2l 
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The  intensity  of  the  beam  is  given  by 


(B2) 


(B3) 


T/-  .  exp[  -  7  V(  C2  +  ( 1  -  5y?)  2)3 
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In  Eqs.  (B2),  (B3),  and  (B4),  "2  =  (x2»yVa2,  l  -  z/1<n2,  and  f  =  f  An2. 

A  collimated  beam  is  obtained  from  the  limit  of  infinite  focal  length: 
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Ic(p,z)  =  f_^-^/(1U2)]  (B7) 

va2(  1  +  £2) 

By  comparing  the  numerically  computed  solutions  for  these  cases  with  the  analytical 
results,  those  portions  of  the  code  that  deal  with  the  linear  parts  of  the  equation  are 
checked,  and  the  linear  stability  and  convergence  of  the  algorithm  can  be  ascertained. 


High-Energy  Beams 

When  the  absorbed  power  in  the  central  portions  of  the  beam  near  the  laser  face  is 
large,  the  nonlinear  term  in  Eq.  (56)  will  be  much  larger  than  the  Laplacian  term.  Such 
a  state  of  affairs  is  local  only,  since  the  opposite  must  hold  in  the  wings  of  the  beam. 
For  this  local  region,  however,  Eq.  (56)  may  be  approximated  by 


From  Eq.  (B8)  and  its  complex  conjugate,  it  may  easily  be  shown  that 


(B8) 
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Hence  the  integral  is  independent  of  the  variable  £,  and  therefore  Eq.  (B8)  may  be  inte¬ 
grated  directly  to  give 
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Equation  (B8)  and  its  solution  (BIO)  prove  useful  in  two  ways.  Regardless  of  how 
good  the  approximation  of  dropping  the  Laplacian  in  any  one  specific  case  may  prove  to 
be,  those  portions  of  the  computer  code  that  deal  specifically  with  the  Laplacian  may  be 
bypassed,  solving  Eq.  (B8)  instead  of  Eq.  (56).  Thus,  a  check  is  provided  on  the  formu¬ 
lation  and  accuracy  of  the  code  in  treating  the  nonlinear  portions  alone.  Also,  for  high 
values  of  the  oscillations  of  the  amplitude  as  a  function  of  ;,  induced  by  the  heating, 
must  be  sampled  at  least  six  times  per  oscillation  to  be  certain  of  accuracy.  Dropping 
the  cooling  portion  of  the  integral  and  using  a  Gaussian  beam,  this  condition  yields  a 
restriction  on  stop  size  of 


.'At,  <  l  (Bll) 

which  was  found,  in  practice,  to  be  necessary  indeed,  as  was  discussed  in  Section  IV. 


Appendix  C 


EFFECTS  OF  THERMAL  CONDUCTION 


In  this  appendix,  the  size  of  the  thermal  conduction  terms,  which  have  not  been  in¬ 
cluded  in  the  analysis,  will  be  estimated,  and  their  neglect  justified.  Thermal  conduc¬ 
tivity,  which  provides  a  mechanism  for  energy  transfer  in  addition  to  the  hydrodynamic 
motion,  changes  the  energy  conservation  Eq.  (7)  to 


dp  , 
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The  effect  of  the  added  terms  is  estimated  using  the  solution  obtained  by  neglecting  them. 
The  pressure  and  its  gradients  are  negligible  in  any  case,  so  only  the  density  terms  are 
considered;  replacing  the  V2  term  by  derivatives  along  the  wind  where  maximum  changes 
occur,  and  pQ  by  a  I,  the  ratio  of  the  conduction  term  to  the  transport  term  then  is 
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The  factor  (l/l)(dl/dx)  determines  a  characteristic  inverse  length  appropriate  to  our 
estimate;  this  length  can  be  estimated  from  the  output  intensity  profiles  at  the  smallest 
cross  section  of  the  beam  and  is  about  1  cm  in  the  most  severe  cases.  Using  the  appro¬ 
priate  thermal  constants  for  air,  and  for  a  wind  speed  of  200  cm/sec,  the  above  ratio  is 
of  the  order  of  10" 3;  when  slewing  is  included,  the  ratio  is  even  smaller.  The  effects  of 
thermal  conduction  in  the  cases  presented  are  negligible  since  this  estimate  provides 
only  an  upper  limit  to  these  effects;  each  case  must,  however,  be  considered  separately. 
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